#include<iostream>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<pthread.h>
#include<time.h>
#include<limits.h>
#define NUM_THREADS 6
#define Nsoma 30

using namespace std;

double *trials;

double S,E,r,sig,T,aux1,aux2,aux3;
unsigned int m;

// gerador de numeros aleatorios

// gera números diferentes a partir do ultimo numero
const unsigned int A = 1664525, C = 1013904223;

unsigned int numale(unsigned int x){
        return x*A+C;
}

// Gerador de numero aleatorio com distribuição gaussiana
double random_number(unsigned int *value){
        double x = 0;
        
        for(int i = 0; i < Nsoma; i++){
                *value = numale(*value);
                x += (double)*value/(double)UINT_MAX;
        }
        
        x -= Nsoma / 2.0;
        x /= sqrt(Nsoma / 12.0);

        return x;
}


double mean(double *v, int M){
        double sum = 0;
        
        for(int i = 0; i < M; i++){
                sum += v[i];
        }
        
        return sum / (float)M;
}

// retorna desvio padrão dos valores de m
double stddev(double *v, float mean, int M){
        double variancia = 0;
        
        for(int i = 0; i < M; i++){
                variancia += pow(v[i] - mean,2);
        }
        
        return sqrt(variancia / (float)M);
}


/* 
	algoritmo de black-scholes
*/
void *black_scholes(void *threadID){

        int taskid, *id_ptr, ini, fim;
        unsigned int num;
        double temp;

        num = time(NULL);
        
        // get task id number
        id_ptr = (int *) threadID;
         taskid = *id_ptr;
        
        // calculate thread trial interval
        ini = taskid * m;
        fim = ini + m;
        
        for(int i = ini; i < fim; i++){
                
                // get trial i
                temp = S * exp(aux1 + aux2 * random_number(&num));
                trials[i] = aux3 * fmax(temp-E,0);
                //printf("trial %d = %f\n",i,trials[i]);
        }
        
        pthread_exit(NULL);
}

int main(){

        double media, desvio, confwidth, confmin, confmax;
        int M,*vetID[NUM_THREADS];
        pthread_t threads[NUM_THREADS];
        
        cin >> S;
		cin >> E;
		cin >> r;
		cin >> sig;
		cin >> T;
		cin >> M;
		
        trials = (double *) malloc (sizeof(double) * M);

        aux1 = (r - (sig * sig / 2)) * T;
        aux2 = sig * sqrt(T);
        aux3 = exp(-r * T);
        
        m = M / NUM_THREADS;

		// executa as threads (N/threads)
        for(int i = 0;i < NUM_THREADS; i++){
                vetID[i] = (int *) malloc (sizeof(int));
                *vetID[i] = i;
        
                pthread_create(&threads[i],NULL,black_scholes,(void *)vetID[i]);
        }
        
        for(int i = 0 ;i < NUM_THREADS;i++)
                pthread_join(threads[i],NULL);

		// calcula os parametros
        media = mean(trials,M);
        cout << media << endl;
        desvio = stddev(trials,media,M);
        confwidth = 1.96 * desvio / sqrt(M);
        confmin = media - confwidth;
        confmax = media + confwidth;

		//imprime dados
        cout << "S\t" << S << endl;
		cout << "E\t" << E << endl;
		cout << "r\t" << r << endl;
		cout << "sigma\t" << sig << endl;
		cout << "T\t" << T << endl;
		cout << "M\t" << M << endl;
	
		cout << "Confidence inteval: (" << confmin << ", " << confmax << ")" << endl;
        free(trials);

        return 0;
}

